| SOCR ≫ | BPAD Website ≫ | BPAD GitHub ≫ |
THIS IS JUST a PLACEHOLDER Template for the BPAD Markdown style….
Let’s examine various strategies to understand and comprehend the complexity of the problem of “optimal” representation of the intricate balance between intervention (LC) suppressing tumor growth and health side-effects of chemotherapy (RP2) in oncology.
The goal is to understand what may be appropriate and simple mechanisms to communicate this balance amongst basic scientists, clinical investigators, and AI algorithms.
In this project, we are using (perhaps over-simplified) exponential models of growth (tumor) and decay (RP2). We’ll simulate realistic scenarios and investigate the univariate patterns and the bivariate relations of LC and RP2. There are lots of parameters and assumptions in this model-based approach. Later we can transition from model-based to data-driven estimates of the parameters to make the optimal surfaces more realistic.
The goal is to answer the following question - "How to make this Data+AI+Forecasting easily communicated (and manipulated) by clinicians and used in patient onco-counseling. As all graphs are really dynamic, the hope is that we may be able to transform classical Age-growth-curves to Optimal-Cancer-treatment-surfaces.
Before proceeding, please review the DSPA R Introduction Chapter.
Suppose for the specific cancer type, tumor cells reproduce at a rate \(r>0\). Then, if \(y_{o}\) is some initial tumor size, the total tumor size (local-control) representing the total number of cells at time \(t\) is given by: \[\underbrace{y\left(t\right)}_{\text{tumor}}=y_{o}e^{rt} .\]
Similarly, we assume that the \(RP_2\) side-effects (pain) increase exponentially with the increase of the chemotherapeutic strength with a rate \(q>0\). Then, if \(z_{o}\) is some initial default pain level, the \(RP_2\) representing the side-effects at time \(t\) is given as a pain value: \[\underbrace{z\left(t\right)}_{{\text{RP}}_2}=z_{o}e^{qt} .\]
Each time a chemotherapeutic agent is introduced, it destroys a fraction \(f\) of the tumor cells that exist at that moment in time. We are trying to develop a model for \(y\left(t\right)\) as a function of time, and the strength of the chemotherapy, which is directly proportional to the size-effects of the treatment applied as \(T\) time separated treatment administrations. We want to explore different cases and explocate the the relations between \(f,T,r\)?
Generally, increasing the chemotherapy dose tends to improve patients tumor LC, however, this highly correlates with increasing risk of RILTs, such as radiation pneumonitis (RP) with grade \(\geq2\) (\(RP_2\)), see this pub PMC6279602.
Build a set of functions:
tumor_run(): the first helper function runs the exponential-growth model each time a chemotherapy administration event occurs, \(y_o \times e^{r \times time},\ r>0\),tumor_redux(): the second function captures the entire duration of the experiment, \(\sum_{k=1}^{Duration/T} {y_k \times e^{r \times time}}\), \(y_{k+1}=y_k\times (1 - f)=y_o\times (1-f)^k\),pain_run(): the third helper function runs the pain exponential-decay model followingeach time a chemotherapy administration event occurs, \(z_o \times e^{q \times time},\ q<0\),painSideEffects_RP2(): the last function tracks the side-effects, pain and RP2, \(\sum_{k=1}^{Duration/T} {z_k \times e^{q \times time}}\), \(z_{k+1}=z_k\times (1 + f)=z_o\times (1+f)^k\).We will examine different scenarios to explore the effects on the clinical outcome (local-control, tumor size) of this model-based simulation:
# function running the exponential model each time a chemotherapy administration event occurs
tumor_run <- function(y0, T, r) {
time <- seq(from = 0, to = T - 1, by = 1)
num_cells <- y0 * exp(r * time)
}
# Aggregate function capturing the entire duration of the experiment
tumor_redux <- function(y0, T, r, f, duration) {
num_runs <- seq(0, duration %/% T - 2, by = 1)
time_tot <- seq(0, duration - 1, by = 1)
cells <- tumor_run(y0, T, r)
for (val in num_runs) {
last_num <- cells[length(cells)]
new_start <- last_num * (1 - f)
new_cells <- tumor_run(new_start, T, r)
cells <- c(cells, new_cells)
}
cells
}
pain_run <- function(z0, T, q) {
time <- seq(from = 0, to = T - 1, by = 1)
pain <- z0 * exp(q * time)
}
# Pain/Side-effects function (returns RP2)
painSideEffects_RP2 <- function(z0, T, q, f, duration) {
num_runs <- seq(0, duration %/% T - 2, by = 1)
time_tot <- seq(0, duration - 1, by = 1)
pain <- tumor_run(z0, T, q)
for (val in num_runs) {
last_pain <- pain[length(pain)]
new_start <- last_pain * (1 + f)
new_pain <- tumor_run(new_start, T, q)
pain <- c(pain, new_pain)
}
pain
}
elapsed = 100
time <- seq(0, elapsed - 1, by = 1)
ctrl_y0 = 10 # starting tumor size (number of cancer cells)
ctrl_T = 10 # treatment-frequency (each T days)
ctrl_r = 0.05 # tumor exponential parameter (growth rate)
ctrl_f = 0.3 # treatment efficacy (fraction of elimination), proportion (1 - f) cells remain after each treatment
treatment <- array(dim=1)
for (i in 1:elapsed) {
treatment[i] <- 0
if (i %% 10 ==0) treatment[i] <- 5
}
# Pain RP2
ctrl_z0 = 5
ctrl_q = -0.05
########## LC
LC <- tumor_redux(ctrl_y0, ctrl_T, ctrl_r, ctrl_f, elapsed)
long_interspersed <- tumor_redux(ctrl_y0, ctrl_T + 10, ctrl_r, ctrl_f, elapsed)
low_reproduction <- tumor_redux(ctrl_y0, ctrl_T, ctrl_r - 0.02, ctrl_f, elapsed)
high_elimination <- tumor_redux(ctrl_y0, ctrl_T, ctrl_r, ctrl_f + 0.3, elapsed)
########### RP2
RP2 <- painSideEffects_RP2(ctrl_z0, ctrl_T, ctrl_q+0.02, ctrl_f, elapsed)
long_interspersedRP2 <- painSideEffects_RP2(ctrl_z0, ctrl_T + 10, ctrl_q+0.02, ctrl_f, elapsed)
heavyCompoundingRP2 <- painSideEffects_RP2(ctrl_z0, ctrl_T, ctrl_q + 0.04, ctrl_f, elapsed)
high_eliminationRP2 <- painSideEffects_RP2(ctrl_z0, ctrl_T, ctrl_q+0.02, ctrl_f - 0.3, elapsed)
df <- data.frame(time, LC, long_interspersed, low_reproduction, high_elimination,
RP2, long_interspersedRP2, heavyCompoundingRP2, high_eliminationRP2)
# fig <- fig %>% add_annotations(text="ctrl_y0 = 10 # starting tumor size (number of cancer cells)\n
# ctrl_T = 10 # treatment-frequency (each T days)\n
# ctrl_r = 0.05 # tumor exponential parameter (growth rate) \n
# ctrl_f = 0.3 # treatment efficacy (fraction of elimination),
# proportion (1 - f) cells remain after each treatment", ax = 50)
# Plot LC Growth Models
library(plotly)
fig <- plot_ly(df, x=~time, width = 1000, name="Time")
fig <- fig %>% add_trace(y = ~LC, name = 'Controled Growth',
name="Number of tumor cells",
mode = 'lines+markers', type="scatter")
fig <- fig %>% add_trace(y = ~long_interspersed, name = 'High T separation (T+10)',
mode = 'lines+markers', type="scatter")
fig <- fig %>% add_trace(y = ~low_reproduction, name = 'Slow Reproduction (r-0.02)',
mode = 'lines+markers', type="scatter")
fig <- fig %>% add_trace(y = ~high_elimination, name = 'High fraction of elimination (f+0.3)',
mode = 'lines+markers', type="scatter")
fig <- fig %>% add_trace(x = ~time, y = ~55*treatment, type = 'bar', # Chemo Treatments
marker = list(color='gray', line = list(color='gray', width=1)),
opacity=0.2, name="ChemoTreat", text = "Chemo Treatments",
hoverinfo = 'Chemo Treatments')
fig <- fig %>% layout(title = "Tumor Growth (Cell Local Control) w.r.t. Different Chemo Interventions vs. Time (T=10, r=0.05, f=0.3)",
scene = list(xaxis = list(title = "Time"),
yaxis = list(title = "Number of tumor cells") ),
legend=list(title=list(text='<b>Models</b>')))
fig# Plot RP2 Pain Models
library(plotly)
# RP2, long_interspersedRP2, low_reproductionRP2, high_eliminationRP2
fig <- plot_ly(df, x=~time, width = 1000, name="Time")
fig <- fig %>% add_trace(y = ~RP2, name = 'Controlled Pain',
name="Pain RP2",
mode = 'lines+markers', type="scatter")
fig <- fig %>% add_trace(y = ~long_interspersedRP2, name = 'High T separation',
mode = 'lines+markers', type="scatter")
fig <- fig %>% add_trace(y = ~heavyCompoundingRP2, name = 'Compounded Pain (q+0.04)',
mode = 'lines+markers', type="scatter")
fig <- fig %>% add_trace(y = ~high_eliminationRP2, name = 'High fraction of elimination (f=-0.03)',
mode = 'lines+markers', type="scatter")
fig <- fig %>% add_trace(x = ~time, y = ~5*treatment, type = 'bar', # Chemo Treatments
marker = list(color='gray', line = list(color='gray', width=1)),
opacity=0.2, name="Chemo Treatment", text = "Chemo Treatments",
hoverinfo = 'Chemo Treatments')
fig <- fig %>% layout(title = "RP2 (Pain) Model w.r.t. Different Chemo Interventions vs. Time (T=10, q=-0.05, f=0.3)",
scene = list(xaxis = list(title = "Time"),
yaxis = list(title = "Number of tumor cells")),
legend=list(title=list(text='<b>Models</b>')))
figThis function, plot_3d(), plots functions that can be described in this parametric form \(y=f(x,t)\). This will allow us to render 3D parametric manifolds as 2D graphs where with the 3rd dimension \(t\) represented as a dynamic variable using a horizontal slider.
\[\underbrace{y\left(t\right)}_{\text{LC}}=y_{o}e^{rt} .\] \[\underbrace{z\left(t\right)}_{{\text{RP}}_2}=z_{o}e^{qt} .\] Therefore, \[RP_2=z_{o}\left (\frac{LC}{y_o}\right )^{\frac{q}{r}} .\]
xLC <- seq(0,275,length=100)
x_RP2 <- seq(0,25,length=100)
# Let’s first embed the objective function in 3D by cine-animating along the time-axis
RP2RP2_objective <- ctrl_z0 * (xLC/ctrl_y0)^(ctrl_q/ctrl_r)
plot_ly(x=~xLC, y=~RP2RP2_objective, mode = 'lines+markers', type="scatter") %>%
layout(title = "RP2=f(LC) (T=10, r=0.05, q=-0.05, y0=10, z0=5, f=0.3)",
scene = list(xaxis = list(title = "LC"), yaxis = list(title = "RP2") ))# plot_3d <- function(f, x, t) {
# x_range <- x
# x_label <- "x"
# slide_range <- t
# frames <- list()
# plt <- plot_ly()
# for (i in 1:length(slide_range)) {
# slide_val <- slide_range[i]
# # Generate the w values for f(x,t)
# w <-(matrix(apply(expand.grid(x_range,slide_val),1,f),length(x_range)))
#
# # A the start, only make the first frame visible
# visible <- i==as.integer(length(slide_range)/2)
#
# # Add this trace (first frame) to the plot
# plt <- add_trace(plt, x=x_range, y=w, mode='lines+markers', type="scatter", visible=visible,
# name=as.character(slide_val), showlegend=FALSE,
# colors = colorRamp(rainbow(8)), opacity=0.5, hoverinfo="none")
#
# # Configure this step in the slider to make this frame visible and none other
# step <- list(args = list('visible', rep(FALSE, length(slide_range))),
# method = 'restyle', label=paste0("t=", round(slide_val,3)))
# step$args[[2]][i] = TRUE
# frames[[i]] = step
# }
#
# # Show the plot + slider focused on the middle plot
# plt %>%
# layout(
# title = paste0("Time-Dynamic 2D Plot LC vs. RP2"),
# scene = list(yaxis=list(title="w=f(x,t)"), xaxis=list(title=x_label)),
# sliders = list(list(active = as.integer(length(slide_range)/2),
# currentvalue = list(prefix = "t:"),
# steps = frames))) %>%
# hide_colorbar()
# }
#
# RP2_objective <- function(LC) {
# return(ctrl_z0 * (LC/ctrl_y0)^(ctrl_q/ctrl_r) ) # see the earlier definitions of the parameters
# }
#
# plot_3d(RP2_objective, xLC, time)# fig <- plot_ly() %>%
# add_trace(data = df, x = ~time, y = ~LC, z = ~RP2, type="mesh3d",
# contour=list(show=TRUE, color="#000", width=15, lwd=10),
# opacity=0.5, hoverinfo="none") %>%
# # trace the boundary of the saddle point surface
# add_trace(x = ~time, y = ~LC, z = ~RP2, type="scatter3d", mode="lines",
# line = list(width = 1, color="red"), name="Optimal Treatment Surface",
# hoverinfo="none")
# fig
fig <- plot_ly() %>%
add_trace(data = df, x = ~RP2, y = ~LC, z = ~time, type="mesh3d",
# colors = colorRamp(rainbow(8)),
opacity=0.7, name="Optim-Treat-Surf", # hoverinfo="none",
contour=list(show=TRUE, color="#000", width=15),
hovertemplate = 'RP2: %{x:0.2f}\nLC: %{y:0.2f}\nTime: %{z:0.0f}') %>%
# trace the boundary of the Optimal Treatment Surface
add_trace(x = ~RP2, y = ~LC, z = ~time, type="scatter3d", mode="lines",
line = list(width = 1, color="red"), name="Opt-Surf-Trace",
hovertemplate = 'RP2: %{x:0.2f}\nLC: %{y:0.2f}\nTime: %{z:0.0f}') %>%
layout(title = "Optimal Treatment Surface", showlegend = FALSE,
scene = list(
xaxis = list(title = "RP2"),
yaxis = list(title = "LC"),
zaxis = list(title = "Time")
))
fig